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We study a statistical mechanics model of a solid. Neighboring atoms are connected by Hookian 
springs. If the energy is larger than a threshold the spring is more likely to fail, while if the energy is 
lower than the threshold the spring is more likely to be alive. The phase diagram and thermodynamic 
quantities, such as free energy, numbers of bonds and clusters, and their fluctuations, are determined 
using renormalization-group and Monte-Carlo techniques. 
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I. INTRODUCTION 

The mechanical properties of sohds, such as the mechanical failure, are topics of considerable interest .^'^Illllll Jn 
this paper we analyze an equilibrium statistical mechanics model^of a solid. In previous calculationd^ve went beyond 
the ideal Hooke law for springs by using the realistic anharmonic energy versus atomic distance developed and 
tested extensively by Ferrante and his collaborators.!^ We found that the phase diagram exhibits universal features 
when the temperature and the stress are appropriately scaled. Those calculations were mean-field in character as 
we assumed that all springs have the same strain. In this paper we evaluate the role of thermal fluctuations by 
using renormalization-group and Monte-Carlo simulations. The model is deflned in Section II. We view the solid as 
a collection of harmonic springs. If the energy of such a spring is larger than a threshold, the spring is likely to 
fail.lSI Assuming that the relaxation times are short compared to the measurement time, we use equilibrium statistical 
mechanics to compute the various thermodynamic quantities. When the elastic energy is not too large, th e part ition 
function for the harmonic "springs" defined on percolation clusters can be m apped into a Potts modeL^^Ull This 
model is quite similar to the annealed Ising model on percolation cluster.^i^liSI Section III we present numerical 
results based on the renormalization-group Migdal-Kadanoff scheme which, as first sho wn by Berkei^^^, provides exact 
solutions of statistical models on hierarchical lattices-'^^ and on small world nets.'^^'i^ Monte-Carlo simulations 
are presented in Section IV. Our concluding remarks are found in Section V. 

II. MODEL 

The energy of a "spring" («, j) is given by the Hooke law: 

m,^-Ec+\{n-r,Y (1) 

where is the displacement vector from the equilibrium position of atom z, measured in units of the lattice constant 
a, Ec is the cohesive energy, and k is elastic constant. If the energy of the spring is larger than the threshold energy 
Eq the "spring" is more likely to fail than to be alive, p is the probability that the "spring" is alive and 1 — p the 
probability that the "spring" breaks. We assume its dependence on energy to be given by the Boltzmann weight 

^^e-'^=z«e-f(--'-^)^ (2) 
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where: K — k/kBT and w — e ''b'^ . We allow for correlations between failing event s by using the Potta^ number 
of states q, which plays the role of a fugacity controlling the number of clusterd^HHIl if g >> 1 there is a tendency 
of forming many small clusters while if g << 1 there is a tendency to form a few large clusters. If q = I springs fail 
independently of one another, i.e. random percolation probleirPSl. The partition function is obtained by summing 
over all possible configurations of bonds arranged on the lattice 

con fig 

C is the number of clusters, including single site clusters, and B is number of live "springs". The restricted partition 
function associated with the elastic energy for a given configuration of bonds (live " springs" ) is 

„config „ - elastic 
^elastic = ^'^'•^ " 

In Eq. ^ the sum is over all live "springs". 

By using the Kasteleyn-Fortuin cxpansioEp^l for Potts model we can rewrite the partition function as 

Z = Tr„Trre^'^ (6) 

The Hamiltonian is 

~ " E ^o) - y<^(f^»: - r^f] (7) 

^ <ij> 

where ui is a Potts spin taking q values. This mapping is a Gaussian approximation valid when, on the right hand 
side of Eq. [7] the elastic energy is small compared to the first energy contribution. The coupling constants Ji and J2 
are related to the original parameters, w and K , as follows: 

Ji = ln(l + it;) (8) 
J2 = (9) 

The free energy per bond is: / = In -^/-^bonds' "^^"^ derivatives of the free energy / with respect to the parameters 
w, K, and q provide respectively the number of live "springs" 6, the elastic energy ^'elastic ^^'^ number of clusters 
c, each normalized by the total number of lattice bonds: 

^elastic = (10) 
" - »-£ (11) 

The derivatives of those densities, ^'elastic ^' ^"^^ ^^^^ respect to the model parameters provide in turn the 
fiuctuations (variances) of those quantities: 

^^llastic - -^'^+Wic (13) 

Ab^ - (14) 
ow 

ic» = .| (15) 
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III. RENORMALIZATION GROUP 



The Migdal-KadanofF recursion equation d^'^l^'^ l for d dimensions are: Zi,j — {TrkZi,kZk.j)^ , where L — 2^d~ 1). We 
assume each atom coordinate varies in the interval (-1/2, 1/2), where the equilibrium lattice constant is 1. After also 
using the Gaussian approximation (small elastic energy) we get: 



w' = [l + U{w,K,q)]^ ~1 (16) 
L 
2 



K'w' = K^[l + U{w,K,q)]^-^U{w,K,q) (17) 



where L = 2'* ^ and 



u('w,K,q) = — , — , (18) 

q^/K/^+V8w.erf{./K/8) 



The recursion equations (16l-(17l represent the Gaussian approximation of the exact solutions for hierarchical lat- 
tices. Since this scheme is realizable, the convexity of the free energy is preservecP^, and thus reasonable expectations, 
such as positivity of energy fluctuations, are fulfilled. The renormalization group flows are governed by the following 
fixed points a,t K ^ (pure Potts model): i. w = (non-percolating live "springs"), ii. w — 8 (percolating network 
of live "springs"), iii. w — Wc (Potts critical point). A stability analysis at the Potts critical point, {K = 0, w = Wc) 
yields the two eigenvalues: i. the thermal eigenvalue Ai (for the direction along the K = axis) is always larger than 
1, meaning the w ^ Wc is a, relevant field; ii. The other eigenvalue A2 is associated with the flow along the w = Wc 
line away from the pure model {K = 0). For d = 2, A2 < 1 for all q. This means that there is a line of points in the 
(w,K) flowing into, and thus is in the same universality class as, the pure Potts critical point (wcO). For d = 3 on 
the other hand, A2 < 1 for q < 109, but A2 > 1 for q > 109. There exists another fixed point at {w*,K*) which has 
both eigenvalues larger than unity for q < 109, and becomes stable in one direction for q > 109. Thus in d = 3, for 
large enough q, the elastic constant K becomes a relevant field changing the universality class of the model from the 
pure Potts criticality to a new one. Potts-elastic. However, in view of the Gaussian approximation used to derive the 
recursion equations, we view this as only an indication of a possible new universalty class that warrants further study. 

The phase diagram for any given q, in the (w,K) plane shows two phases: I. solid with a percolating network of 
live "springs", II. crumbling solid with mostly "failed" springs. The two phases are separated by a critical line in the 
universality class of the g-state Potts model (for d — 2 for all q, for d = 3 for q < 109). 

Note that when increasing the elastic constant K, one needs a higher w to establish the solid phase. This is due to 
the fact that increasing the elastic energy increases the probability for the "spring" to fail. 

The free energy / — In Z/Nb , where Nb is the number of lattice edges, is: 



- {2Ly 

where 



(19) 



C = Hq + w^-erf{^-)]^ (20) 

and L = 2'^~^ . Using the free energy we can compute the number of live "springs" 6, the number of clusters c, the 
elastic energy ^'elastic their fluctuations (variances). Each of those quantities is scaled by the total number 
of lattice edges Nb- In Figure 3 we show the elastic energy variation with K and w for two different q values. As 
expected the elastic energy increases monotonically with K and w starting at zero at K = and at w = 0. 
The number of live "springs" increases with w and decreases with K, as shown in Figure 4. 

We also estimate the squared mean elongation (in units of lattice constant a) of live " sprins" by using the number 
of live springs, 6, and the elastic energy: 

Ar' ^ j_ -^elastic (21) 
K b 

One can use the classical Lindemann model of meltin^^, Ar^ = 0.01, to estimate the melting temperature of our 
model solid. 
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FIG. 1: Thermal and Elastic Eigenvalues in d = 3. 



FIG. 2: Phase diagram for q = 10 and d = 2, in the plane {w, K). 



In the limit w — Q, the number of clusters is equal to the number of sites. Hence c approaches the inverse of the 
coordination number, which for the diamond hierarchical lattice, correponding to the Migdal-Kadanoff scheme for 
d = 2, i^^. c = 2/3, consistent with Figure 6. 

In Figure 7 we show the elastic energy fluctuations and the number of live "springs" fluctuations as functions of 
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FIG. 3: Elastic energy vs K and w. Left column q — 10; Right column q — 1. 




FIG. 4: Number of live springs vs. w and K, respectively. Left column q = 10; Right column q — 1. 



K and w, for q = 10 and g— 1 respectively. Since the exponent a is positive for q — \Q and negative for q — I, a, 
divergence is apparent in the q = 10 critical point K ^ 2, w ~ 7.3. 

The number of clusters c increases monotonically with the conjugated fugacity q, starting at c = at q = 0, as 
shown in Figure 8. The fluctuations in c exhibit a divergence at the critical point <; = 10, K— 2, w — 7.3 where the 
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FIG. 5: Square of live "spring" elongation vs. w and K, respectively. Left column q = 10; Right column q = 1. 




FIG. 6: Number of clusters vs. w and K, respectively. Left column q = 10; Right column q — 1. 



critical exponent a is positive. 

In Figure 9 we show fluctuations in elastic energy, in number of bonds, and number of clusters against w, K, and 
q respectively, ford = 3 and q = 20. One can notice the lack of symmetry in the divergence which is a characteristic 
of 3d criticalit} *^^ l ^^ l By contrast in 2d the divergences are symmetric (see Figure 7). This is related to the duality 
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FIG. 7: Fluctuations (variances) in elastic energy and number of live "springs" vs. K and w, respectively. Left column q = 
10; Right column q = 1. 




FIG. 8: Number of clusters and their fluctuations (variances) vs. q. Left column w — 7.3, K — 2; Right column w — 1.857, 
K ^2. 



transformatiorPSEZl The critical point exponent at g = 20, w = 3.9, = 2 is a = 0.015. 
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FIG. 9: Fluctuations of elastic energy, number of live "springs", and number of clusters vs. k, w, and q, respectively, for d = 
3. 



IV. MONTE CARLO SIMULATIONS 



In this section, we use for Monte Carlo simulation the following Hamiltonian taken from Eq. (7) 

H=-h '^j) + t E - ^i)^ (22) 

<ij> <ij> 

where Ii and I2 are renormalized parameters. Of course, one has Ji — Ii/ksT and J2 = 12/kBT. 

In the case of two dimensions d = 2, we consider a square lattice of size N x N where N = 40, 60, 80, 100. Each 
lattice site is occupied by a g— state Potts spin. We use periodic boundary conditions. Our purpose here is to locate 
the phase transition point and establish the phase diagram. To this end, a simple heat-bath Metropolis algorithm is 
sufficient. The determination of the order of the phase transition and the calculation of the critical exponents in the 
second-order phase transition region need more sophisticated Monte Carlo methods such as histogram techniques.'^ 
These are left for a future study. 

The simulation is carried out as follows. For each set of (/i, I2), we equilibrate the system at a given temperature T 
during 10^ Monte Carlo sweeps (MCS) per spin before averaging physical quantities over the next 10® MCS. In each 
sweep, both the spin value and the spin position are updated according to the Metropolis criterion. The calculated 
physical quantities are the internal energy per spin E, the specific heat Cy per spin, the Potts order parameter Q and 
the susceptibility per spin x- For a g-state Potts model, Q is defined as 



Q ^ gmax(Qi,Q2,---,Qq) - 1 n^^. 
q-l 

where Qi = ^ (z = 1, g), being the number of sites having qi. 

Let us show first in Fig. [TOjand Fig. [TTJthe energy, the specific heat, the order parameter and the susceptibility in 
the case where q — 2, Ii ^ 1 and I2 = 0.5. 

These figures show a phase transition at Tc — 1.130 ± 0.005. Note that the size effects for N = 40, 60, 80 and 100 
are not significant and are included in the error estimation. Simulations have been carried out also for the following 
sets (/i = l,l2 — 0.2), (/i = l,/2 = 0.8) and (/i = l,/2 = 1). The results show that the transition temperature Tc 
does not change significantly with this range of 12- Tc depends only on the main Ii term. 
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FIG. 10: Energy per spin E (left) and specific heat per spin Cv (right) vs temperature T for q = 2^ I2 = 0.5 with = 60 and 
h = 1. 
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FIG. 11: Order parameter Q (left) and susceptibility x (right) vs temperature T for q = 2, I2 ~ 0.5 with TV = 60 and 7i = 1. 



To compare with the results from the renormahzation group calculation of the previous section, we have to use 
Ji = Ii/ksT, J2 = h/ksT and Eqs. (8]|9l to convert Ji and J2 into K and w. One has 



w = exp(/3/i) — 1 
h 1 



K = 



keT 1 - exp(-/3/i) 



(24) 
(25) 



where (3 = l/fcsT. 

Figure [12] shows the internal energy E as functions of w and K for Ii — I and I2 — 0.5. The transition is found at 
w* — 1.421 ± 0.002 and K* = 0.753 ± 0.002. We observe that only K* varies with I2, not w* as expected from Eqs. 



(24l-(25l. We have K* = 0.301,1.205 and 1.506 for I2 = 0.2,0.8 and 1, respectively. 



V. 



CONCLUSIONS 



We have studied a model of a solid made of springs that are live and harmonic or failed. The springs can fail with a 
probability that increases with the energy. Our renormalization-group analysis suggests that elastic perturbations on 
the Potts-percolation model are irrelevant for all q in two dimensions, and for small enough q in three dimensions. The 
renormalization-group predictions must be viewed as only indicative, in view of the following known limitations. The 
simple Migdal-KadanofF renormahzation group fails to predict correctly the first-order transitions for d — 2, q > A, 
and for d = 3, g > 2. Furthermore our recursion equations are valid for small elastic energy. Monte Carlo simulations 
are needed to further our understanding of the model. Monte Carlo results in d = 2 show that the phase transition 
does not depend on I2 both on the value of the transition temperature and on the transition order. This is physically 
in agreement with the fact that in two dimensions the melting does not take place at finite temperature. At least in 
d = 2, one can say that the phase transition is solely due to the first term of the Hamiltonian (22 1. We note in passing 
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FIG. 12: Internal energy E versus w (left) and versus K (right) for q = 2, I2 = 0.5 with = 60 and Ii — 1. 



that in the Ising case, the question of whether the elastic intera ction affects or not the Ising universahty class has 
been investigated.'^l^ No definite conclusions have been re ached Jt would be therefore interesting in the future 
to perform Monte Carlo simulations for d = 3 and for large q to investigate the effect of elastic interaction on the 
criticality and on cross-over from second to first-order. 
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